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This paper describes a novel approacli based on "proportional 
imputation" when identical units produced in a batch have random 
but independent installation and failure times. The current problem 
is motivated by a real life industrial production-delivery supply chain 
where identical units are shipped after production to a third party 
warehouse and then sold at a future date for possible installation. Due 
to practical limitations, at any given time point, the exact installa- 
tion as well as the failure times are known for only those units which 
have failed within that time frame after the installation. Hence, in- 
house reliability engineers are presented with a very limited, as well as 
partial, data to estimate different model parameters related to instal- 
lation and failure distributions. In reality, other units in the batch are 
generally not utilized due to lack of proper statistical methodology, 
leading to gross misspecification. In this paper we have introduced a 
likelihood based parametric and computationally efficient solution to 
overcome this problem. 

1. Introduction: Background of the problem. After the production pro- 
cess, consumer goods are often distributed through muhi-step channels, giv- 
ing rise to the term "production-delivery" supply chain. An exception to this 
practice is "just-in-time" manufacturing where a product is assembled and 
shipped directly only upon the request of a customer, which is quite popular 
in the personal computer industry. However, for most consumer products, 
items produced by a company are not shipped directly to the final cus- 
tomer. The traditional route for any large scale industrial operation is to 
ship the manufactured products to a warehouse. The warehouses are often 
maintained by third party retailer /shops, from where the products are sold 
and installed at a future date to the final customer. Due to geographic as 
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well as company-retailer relationship, once the batch is shipped, it is often 
unknown to the producing company whether a specific unit is working or 
is still not installed, until and unless the unit stops working and the final 
customer claims a warranty at a future date. At that point in time the 
data on the failed unit becomes "complete" in a sense that we know exactly 
its installation as well as failure time. For all other units it is not known 
(hence "partial" information only) whether they are working or are not at 
all installed. The above setup is quite common in practice in many indus- 
trial supply chains, giving rise to a situation where in-house engineers face a 
dilemma regarding the optimal usage of available information. The untimely 
failure of a unit is always costly to the producer from the warranty perspec- 
tive [Abernethy (1996)]. Also, after infant mortality, reliability assessment 
and future lifetime prediction at an early stage of the product lifespan is 
advantageous for appropriate customer satisfaction issues. 

Reliability estimation requires knowledge of the population at risk and 
the reliability of each unit of the population. The major objective is always 
to acquire timely information of interest on failure modes. However, in the 
presence of both "complete" and "partial" information, current practice is 
to estimate relevant reliability information by using those units which have 
completed their life cycle (i.e., "complete" portion only), while not utiliz- 
ing the "partial" information [Abernethy (1996); Kececioglu (1993)]. The 
primary reason for this is the absence of any established methodology for 
dealing with the current situation. This clearly makes the inferential proce- 
dure suboptimal. In this article we adopt a proportional imputation based 
approach to yield a practical solution to the situation described above. The 
thrust of this paper is the estimation of the unknown parameters under the 
assumption that we know the actual parametric distribution of installation 
as well as failure time. The more general problem of unknown distributional 
form for either installation or failure time (or both) is not considered here 
and is left for future work. 

The rest of the article is organized as follows. In the first three sections 
we present notation and a theoretical justification of the proposed method- 
ology. Section 5 presents the algorithm for proportional imputation. The 
connection between the exact likelihood based approach and our proposed 
algorithm is described in Section 6. Section 7 describes the simulation per- 
formance of our algorithm. We also include the analysis of industrial furnace 
data in Section 8. We conclude the article with some discussion. 

2. Notation and mathematical setting. The problem of interest is moti- 
vated from a large industrial company producing residential furnace com- 
ponents. The units are produced and shipped within the continental USA 
via multiple channels. However, the general description of the problem and 
our solution is neither dependent on a specific company nor confined to a 
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specific commodity. Rather, our proposed solution will have a broader ap- 
plication since the setup is common to many production delivery supply 
chains. Consider a setup in which N identical units are produced in a batch, 
which are then shipped to a warehouse. These units will be installed only 
after being purchased by the customer at some future date. We assume 
there exists no substantial time lag between purchase and actual installa- 
tion of unit /units. Purchase and installation will be considered as the event 
of interest, and the time in which this transpires will be referred to as the 
"installation time." Consider a fixed end of study time Tq. The general data 
description at hand is rather simple. For a particular unit we either know 
both the installation and failure times or know nothing at all. In fact, for 
many units at time Tq, their current status will be unknown due to the 
fact that they have not yet failed either due to noninstallation or are still 
in working condition. Let X ™d T denote the contin- 

uous random variables corresponding to installation time and failure time 
and which are assumed to be independent of each other. In this paper we 
assume that Fx{-) and -Ft(') are completely specified but with unknown 
parameters. We denote the random set i7 = {z S {1, 2, . . . , N} : Xi + Ti < Tq} 
to be the set of indices of the completely observed units. Let C denote the 
cardinality oiQ,:C =\Q\ = "^j^i I{Xi+Ti <Tq}. Following standard results 
in survival/reliability analysis, the complete likelihood for the above setup 
is 

N 



L{Fx,Ft) = ll[fx,T{x^,ti)I{x, +U< To}Y' [P{X + T> To}]'-^^ 
1=1 

(2.1) 

- w .To -^N-C 

Ylfx,T{^uti)j^Sx{To) + J ST{n-x)dFx{^ 



o^<. \_\_Jx.T(Xi,ti]-f<, bx[ln]+ I bT[±n-x]atx[x] 

where Tj is an indicator of whether the ith unit is observed or not for i = 
1,2, ... ,A^. The above likelihood is difficult to maximize numerically except 
for the very restrictive case when X and T are independent and identically 
distributed (i.i.d.) according to an exponential distribution. For the other 
popular reliability distributions (e.g., Weibull, Gamma), the above likelihood 
is difficult to maximize due to excessive fiatness, especially when C <C A^. In 
the furnace data described in Section 8 and also in other simulation studies, 
the ratio is on average 40% or below. With only this much data the 
above likelihood essentially becomes very flat and brute force optimization 
often produces unstable estimates with large variances. For more details on 
this see the simulation studies in Section 7. Next we provide a proportional 
imputation scheme that has close connection with the above likelihood, yet 
it employs a search strategy parallel to Monte-Carlo-based approaches which 
is computationally faster and produces stable estimates. 
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2.1. Standard practice and an alternative formulation. For notational 
simplicity and without loss of generality, we assume that the first C units 
are observed or, in other words, we have complete information for 
Notably, the manufacturer knows nothing about a unit under two circum- 
stances. First, if X > To, that is, the unit is not being installed until time Tq 
and denoted as event B. Second, X <Tq but T > Tq — X, that is, the unit is 
installed but still in operation and denoted as event D. Since exact likelihood 
is difficult to use, traditional practice is of two forms [Abernethy (1996); Ke- 
cecioglu (1993)]. The most simplistic approach is to think that only C units 
are produced. Since we will have complete information for all of them, we 
may use standard theory to estimate model parameters corresponding to X 
and T under specific distributional choices. The other practice is to think 
that we have C units not from the full distribution but rather from the trun- 
cated distribution of both X and T (i.e., observed if X < Tq and T <Tq). 
Then under some specific distributional assumptions (popular choices are 
Exponential, Weibull, etc.) the MLE or rank egression based approaches are 
used for parameter estimation [Wang (2004); Johnson (1964); Michael and 
Schucany (1986)]. Both of these approaches will produces erroneous esti- 
mates for the setup considered. The situation will be much simpler if it is 
also known for a specific "noninformative" unit whether it is under the event 
B or D. This knowledge, if available, will enable us to render the case as 
Type-1 right censoring at Tq either on X (under B) or on T (under D) and 
then follow the usual theory of estimation with censored data [Meeker and 
Escobar (1998); Klein and Moeschberger (2005)]. Unfortunately, practical 
considerations suggest that even this information will not be available un- 
der most producer-retailer setups resulting in "ambiguous" censoring. This 
is unavoidable unless the producer company has an agreement with the re- 
tailer to get in-time unit specific sales information. This involves monetary 
implications and often short-term cost cutting actions get higher priority. 

In this article we took an alternative route to impute the installation 
time (X) for those units under D, that is, installed but not failed. Note 
that if we know or can successfully impute the installation time and assume 
that the unit is still working, this essentially means the failure time is be- 
ing censored. This enables us to use standard methodology to estimate the 
model parameters [see Meeker and Escobar (1998); Klein and Moeschberger 
(2005)]. However, the crucial question is not only how to impute the unob- 
served installation time, but also how many units are needed to be imputed. 
Next we present the theory of an interesting computational approach to 
achieve this task based on a proportional sampling imputation scheme. 

3. How many to sample and where to sample from? In the parametric 
setup we generally assume some distributional form for X and T, Weibull 
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and Exponential being the most popular choice to reliability engineers [Aber- 
nethy (1996)]. Our present methodology is general in the sense that it does 
not depend on any specific distributional choice for both X and T. Note that 
for C complete units we have samples from three conditional distributions, 
namely: 

1. x\X + T<To; 

2. t\X + T<To; 

3. x + t\X + T <To. 

It is not difficult to formalize an estimation procedure if we have samples 
from {x\X < To}. However, the identity 

fx(x\X < n)FT{To - x)Fx(To) 



fx{x\X + T<T, 



implies 



fxix\X<To) 



Ft+x{To) 
fx{x\X + T<To)FT+x{To) 



Ft{To-x)Fx{To) 
(3.1) cxfx{x\X + T<To)F^\To-x) 

oc fxix\X + T< To){l - St{To - x)}-\ 

Remark. Note that the number of samples (if available) from < 
To} will be larger than that from {x\X + T < Tq}. Hence, we have the iden- 
tity, # samples {x\X < To} — # samples {x\X + T < Tq} = ^ samples {x\X < 
To n T > To — X}. We will try to impute this difference (or unobserved in- 
stallations) via proportional sampling. 

The above calculation shows why the assumption that the samples are 
from right truncated and independent distributions is not valid. Even though 
X and T are assumed to be independent, the very nature of the "installation- 
failure" setup will make them intrinsically dependent. Hence, it will be wrong 
to carry out separate estimation of the parameters of the distributions of X 
and T under the truncation assumption, as in reality we do not have sam- 
ples from {x\X < Tq} and {t\T < Tq}. Next we have exploited this mutual 
dependence of X and T via a sampling and imputation based approach. 

3.1. Proportional imputation scheme. To estimate the number of impu- 
tations necessary, let us denote the random variable V = X^^Li where 

1, if jth unit is installed on or before To, 
0, otherwise. 



Hence, P[Vj = I] = P[X < Tq] = Fx{Tq) and Vj ~ Bernouni(Fx(ro)). 
Under the assumption that units are identical and independent, V ~ Binomial (A^, 
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Fig. 1. Schematic diagram of C observed installations. 



Fx(Tq)). Hence, E[V] = NFx{Tq) and since C units are already observed, 
we need to impute for NFx{To) — C units. Of course, NFx{Tq) — C need 
not be an integer and so we round it up to produce a sensible estimate. 
We use [•] notation to denote this rounding procedure. All these make sense 
provided we know the parameters in Fx{-)^ but, in fact, the main purpose 
of this paper is to estimate those parameters. However, for the time being 
let us assume that some crude estimates of these parameters are available. 
We will describe exactly how to get such accurate estimates in Section 5. 

Without loss of generality, we assume C units are ordered in the sense 
that Xi < Xj+i for i = 1, . . . ,C — 1. The observed installations are depicted 
in Figure 1. These installations produce a natural C + 1 partitioning of the 
study interval, that is, [0,To]. Due to the continuous distributional choice 
for X, we consider the case with no ties. However, we remark that the case 
with ties can be handled with minor modifications. The probability of a 
unit being installed in the interval [x^, Xfc+i] is given by P[xk < X < Xk+i] = 
Fx{xk+i) — Fx{xk)- An installed unit will remain unobserved if it does not 
fail by Tq. So the conditional probability of remaining unobserved is given 
by 

, , , , , r'+'ST{T,-x)fx{x)dx 

(3.2) P[T > To -XI.. < X ^ V.(..,0-^xfa) • 

Next we present a theorem for the above conditional probability if the in- 
terval becomes narrower, that is, x.+ilrr.. 



Theorem 3.1. lim,,^,^,, F^Zl^Fl^tl) ^ = ^^^^^ " ^^^^^^^^ 



ST{To-x)Jx{x)dx 

fx{xk+i) /O. 

Proof. This follows by application of I'Hospital's rule. □ 



Remark. This indicates that if Xk+i ixk-, then the probability of sur- 
vival (i.e., remaining unobserved) for a unit installed exactly at Xk will be 
SriTQ-Xk). 

Now using equation (3.2), the joint probability of a unit being installed 
in and then remaining unobserved is 

(3.3) P[{xk <X< Xk+i) n (T > To - X)] = St{To - x)fx{x) dx. 

Jxk 
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Due to the nonincreasing property of the survival function, it is easy to see 
that 

SriTo-Xk) / fx{x)dx< I St{To - x)fx{x) dx 

(3.4) 



< ^^(ro - Xk+i) / fx{x) dx. 

We would like to use the above inequality to approximate equation (3.3) via 

4+1 = P[{xk <X< Xk+i) n (T > To - X)] 

^ ^"^"^"""^^^Y" ^'^°""^"^^ [Fx(x.-,,) - Fxix,)]. 

Remark. Note if Tq J, but [xfc, x^+i] remains fixed with Xk+i < To, then 
Ifc+i t due to the monotone decreasing property of the survival function. 
Conversely, if Tq '[, then Ik+i i- The approximation for Ik+i given in equation 
(3.5) works very well provided the observed installation times are not very 
sparse over [0,To]. Next, we present a theorem characterizing unobserved 
installation times over different regions. 

Theorem 3.2. Let Xk € {xk-i,Xk+i). Then P[T > Tq - X\xk~i <X < 
Xk] < P[T > To - X\xk <X< Xk+i] . 

The proof is provided in the Appendix. Theorem 3.2 implies that the 
probability of remaining unobserved increases as the installation time gets 
closer to the end of study time To. Equation (3.5) characterizes the probabil- 
ity of a single unit being installed in [x^, but remains unobserved until 
Tq. Note that we have C + 1 such intervals in [0, To]. Hence, the expected 
number of unobserved installations in [xk,Xk+i] is 

_ {NFx{TQ)-C}Ik+i 
"fe+i c J ' 

Z^j=Qj^j + l 

with the identity Yl^=Q'^k+i = XFx{Tq) — C. 

Lemma 3.1. X:fc=o4+i = ELi ^4^['5t(To-x,_i)-5t(To-x,+i)] + 
Fx (To) ^^^t{To xc) ^ y^^gj^g xq = and xc+i = Tq. 

Proof. Note that h+i = ^HTo-x.,)+5t(To-x,+0 [^^(^^^^) _ Fx{xk)]. 
Hence, 

E^^+i = E ^"^^^^ " "'^ " "'^'^ iPxixk-,,) - Fxixk)] 

k=0 k=0 
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[Fx{xi) - Fx{xo)] 



St{To - xo) + St{To - x^) 



2 



+ [Fx{x2)-Fx{xi)] 



St{To - xi) + St{To - X2) 



+ [Fx{xc+i)-Fx{xc)] 



St{To - xc) + St{To - xc+i) 



2 



After cancelling successive terms and setting 5^(0) = 1, we complete the 
proof. □ 

Note that even if the distributional forms for X and T are known, ak+i 
will still not be available if we do not know the parameters of Fx{-) and 
Ft{-)- In Section 5 we will propose a general iterative approach for estimating 
these parameters which in turn will yield the estimate a^+i for k = 0, . . . ,C . 
In practice, we use [Sfc+i] for obvious reasons. We would like to put forward 
a sampling based approach to impute these unobserved installation times 
in Section 5. We denote the random set F = {i € {1, 2, . . . , N} : {Xi < Tq) n 
{Xi + Ti> To)} with \T\ = Y.'i:^ o[o^fc+i] being the number of imputed samples 
of X. In this situation, by combining the observed and imputed samples we 
have the case of type-1 right censoring for the installation time X. The 
likelihood for X is then given by 



which we need to maximize with respect to the parameters to obtain the 
ML estimates. 

4. Characterization of failure time. So far our effort was to characterize 
the expected number of unobserved installation times in different partitions 
of [0, Tq] . Once this is known, we want to impute these installation times in an 
iterative fashion (see Section 5). For the time being, if we assume the imputed 
samples represent the actual unobserved installation times, it presents the 
case of random right censoring for T. This is explained in Figure 2. The 
left-hand diagram in Figure 2 represents the possible scenarios with both 
installation and failure times. In the right-hand diagram of Figure 2 we plot 
the time to failure for each unit, taking installation time as the starting 
point. For the imputed installation time (i.e., unobserved due to the fact 
that the unit is still working) what we really get is Tq — X or the random 
censoring time. Hence, the observed variable is T* = minjT, Tq — X}. Note 
that X and T are assumed to be independent and so are T and Tq — X . Let 5 
indicate whether T* is censored {5 = 0) or it is a real failure {6 = 1). For the 



(3.6) 
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Unitl 
Unit2 
Units 
Unit4 
Uniti 
Units 
Unit? 



Installation Time + Time tO' 
failure 



Time to failure 



Fig. 2. Schematic diagram (on left) until observation time To with N — 7 and C = 3. 
A "•" indicates an installation and a indicates a failure. A solid line indicates an 
observed unit (i.e., X + T <To). A dashed line indicates an unobserved unit [i.e., either 
{X > To} or {{X < To) n (T > To - X)}]. Note that units 1, i and 6 are installed but still 
working, while unit 2 is not installed at all. The diagram at the right indicates the time to 
failure only starting from the installation time for each unit (starting from •, at the left). 
Unit 2 does not appear on the right diagram as it has not been installed yet, while units 1, 
4 and 6 are censored for T. 



current situation we have C real failures and [NFx{Tq) — C] censored times, 
while [A^(l — FxiTo))] units do not contribute to the estimation process as 
they provide no information related to failure. The data from n = [N Fx{Tq)\ 
units consists of the pair Since we are interested in inference about 

the parameters of Ft{-), the likelihood function for the same is given by 

n 

(4.1) LT = ll[fTmHST{t*)]'-'^. 

1=1 

5. Iterative algorithm. All our earlier calculations are solely for the pur- 
pose of parameter estimation in the distributions of X and T. The key quan- 
tity of the whole discussion is a^+i (see Section 3.1), which represents the 
number of unobserved installation times in However, the estima- 

tion of ctk+i requires knowledge of the parameters in the distributions of 
X and T. We have assumed so far that the distributions of X and T are 
known; however, the parameters are actually unknown. Hence, an iterative 
procedure is proposed. 

Begin procedure 

Step 0. Find initial parameter estimates of Fx{-) and Ft{-) assuming that 
they are coming from a truncated distribution (< Tq) for which we 
have complete knowledge (e.g., Weibull, Exponential, etc.). 

Step 1. Using the current value of the distribution parameters, find Q?fc+i 
for k = 0, . . . ,C . Note that it is quite possible to have Sfc+i not as 
an integer, say, Sfc+i = int(Sfc+i) + frac(Sfc+i) = Uk+i + Vk+i- 
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Step 2. Draw Uk+i samples from the interval [x^, x^+i] of the distribution 

Fx{-) using current values of the distribution parameters. 
Step 3. First, draw a sample from a Bernoulli(Vfc_|_i). If it is equal to one, 

draw another sample as in step 2, otherwise skip to the next step. 

Hence, the total number of imputed samples is either Uk+i or Uk+i + 

1. 

Step 4. Re-estimate the parameters of X using both imputed and observed 

(C) samples via MLE under right censoring using equation (3.6). 
Step 5. Re-estimate the parameters of T by using both observed (C) and 

censored samples via equation (4.1). The random censoring value 

for any imputed sample is Tq — Ximputed • 
Step 6. Return to step 1 until an acceptable convergence tolerance level is 

reached on the parameter estimates. 

End procedure 

Note that the conventional approach stops at "Step 0" without any fur- 
ther iteration, so we are simply using that as the initial guess. Details for 
obtaining the MLE for some of the truncated distributions (e.g., Exponential 
and Weibull) are described in the Appendix. Though this algorithm assumes 
that the parametric form of X and T are known, it does not depend upon 
any specific distributional choice. Under the assumption that the specific dis- 
tributional choices of Fx{-) and Ft{-) are correct, the speed of convergence 
depends upon the actual observed sample size (C) and end of study time 
(Tq). If C is too small, it will require many imputations (as [NFx{Tq) — C] 
is big). Similarly, if Tq is too small thus representing an early study termi- 
nation, it will force C to be quite small. Both of these cases represent very 
little available information. This generally results in large sampling variance 
with high fluctuations in the iterations resulting in nonconvergence. 

6. Connection with the exact likelihood. Note that our main goal is to 
estimate parameters in the distribution of X and T and typically a likeli- 
hood is a function of those parameters. As noted earlier in Section 3, though 
X and T are assumed to be independent, the nature of ambiguous censoring 
make their joint distribution dependent, where the functional component 
related to respective parameters are nonseparable. As a consequence, max- 
imum likelihood estimation requires joint maximization for all parameters 
over the exact likelihood function given in equation (2.1), which is compu- 
tationally prohibitive. Thus, a major point in this article is the separation 
of the X and T distributions via equations (3.6) and (4.1). A pertinent 
question is the theoretical justification of the above in light of the exact 
likelihood. Note that P{X + T > Tq} = Sx{To) + f^" St{To - x) dFx{x). In 
case there is an oracle which supplies us information about the N — C un- 
observed units, that is, whether {X > Tq} or {X < Tq} n {T > Tq — X}, the 
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above expression simplifies considerably. Suppose that out of those N — C 
units we know that |r| (~ [NFx{Tq) — C]) units are installed (with re- 
ported installation times) but have not yet failed by Tq; then for those 
units, P{X + T > To} = fx{x)ST{To - x). For the remaining iV - |Q| - |r| 
(~ [N{\ — Fx{Tq))]) no information is available, as they are not installed. 
Hence, we get type-1 right censoring on X at Tq, implying P{X + T > Tq} = 
Sx{Tq). The likelihood contribution from the imputed and unobserved units 
is {W-^Y fx{xj)ST{TQ — Xj)}Sx{TQ)^^\^\~\^y Under the above setup, the 
complete likelihood for all observed and imputed samples becomes 

L{Fx,Ft) 

(6.1) oc|n/x(x.)/r(io||n/^(^^)'5T(ro-x,)|5x(ro)^"l^l^l^l 

oc|5x(To)^-l^l-l^l n fx{x^)\{'Y{fT{U)\{ST{To-x,)\. 

This is what corresponds to equations (3.6) and (4.1). 

7. Simulation studies. Next we present some simulation studies with dif- 
ferent choices of reliability distributions to demonstrate the efficacy of the 
proposed approach. In particular, we consider exponential and Weibull dis- 
tributions for both X and T with different values of Tq. To explain the con- 
vergence criteria let us assume is a parameter (in either X or T) that needs 
to be estimated. We stop the iteration when \ti±£-J!±\ < where i denotes 
the iteration number, p is a prespecified positive integer constant and e is a 
prespecified small value chosen by the end user. For multi-parameter cases 
this needs to be satisfied for every parameter. Alternatively, in the spirit of 
the Monte-Carlo-based approach, we may run a fixed but large number of 
iterations and discard the first few iterations as nonstabilized (or "burn-in") 
values and keep all the remaining to report the estimated empirical mean 
and standard deviation. We took the second approach as we found that 
convergence is very fast even for e = 0.0005, except for the situation when 
^ < 20%. In every situation we also report the exact stopping time if we 
choose to use the first stopping criterion (i.e., stop if |/fi±£-Jfi| < s). We also 
report the exact runtime in every simulation using R code on a Windows- 
XP-based machine until convergence. We hope this should give the reader a 
comprehensive idea about the run time efficacy of our approach. The com- 
puter code used for the simulation is available as a supplementary material 
[Ghosh (2009)]. 

Table 1 represents the simulation results for different choices of distribu- 
tions for X and T. We choose = 200 for all experiments. We run the iter- 
ation 1000 times for each model, of which we discard the first 100 as burn-in 
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X ~ Exp(A = 0.5) 


6 


108 


84 


A = 0.69 


A = 0.48, 5a = 0.018 


78 


32 


133 


T~Exp(5 = 0.2) 








(5 = 0.4 


?= 0.22, 5a = 0.01 








X ~ Exp(A = 0.5) 


4 


66 
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A = 0.81 


A = 0.43, 5a = 0.025 


98 


65 
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T ~ Exp(<5 = 0.2) 








5 = 0.56 


5 = 0.23, 5a = 0.013 








X ~ Exp(A = 0.4) 
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A = 0.44, 5 = 0.013 
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108 


T ~ Exp(5 = 0.7) 








5 = 0.78 


5 = 0.7, 5a = 0.03 








X - Exp(A = 0.4) 


4 


124 


43 


A = 0.67 


A = 0.41, 5 = 0.03 


36 


212 


155 


T ~ Exp((5 = 0.7) 








5 = 1.07 


5 = 0.75, 5a = 0.06 








X ~ Exp(A = 0.7) 


6 


111 


84 


A = 1.05 


A = 0.66, 5a = 0.03 


88 


66 


123 


T ~ Weibull(/3 = 2, 6» = 5) 








/3= 1.21 


/3 = 2.03, 5;3 = 0.04 








P = Shape, 9 = Scale 








6' = 1.71£ + 03 


= 5.04, 5e=0.14 








~ Weibull(/3 = 1.5,e' = 4) 


6 


107 


47 


/J = 2.51 
= 5.15 


P = 1.51, 5;3 = 0.04 
= 3.38, 5s =0.14 


35 


45 


115 


T ~ Exp(A = 0.5) 








A = 0.74 


A = 0.44, 5a = 0.033 
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X estimate Plot T estimate plot 




400 600 
Iteration 



(c) 



(d) 



Fig. 3. Plot of the maximum likelihood estimate over different iterations for nonidentical 
exponential cases. Plots at the top are for A = 0.5, 5 — 0.2 and at the bottom are for A = 0.4, 
S = 0.7. The "■ ■ ■ " (dashed line) indicates the true value of the parameter m each case. 



values. The reported parameter estimates and standard deviations are based 
on the remaining 900 iterations. We also report the convergence iteration 
number, which, for the multi-parameter case, represents the maximum of all 
iterations taken by individual parameters to satisfy | ^' I < As we can 
see from Table 1, convergence is achieved quickly. For parameter estimation 
we used the maximum likelihood approach which is described briefly in the 
Appendix section. Again for other nontrivial distributions with complicated 
MLE, the method of moments or rank regression based approaches [Johnson 
(1964)] could be used. In each model, following standard practice, we obtain 
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X estimate Plot T estimate plot 




Iteration Iteration 

(a) (b) 

Fig. 4. Plot of the maximum likelihood estimate over different iterations for the i.i.d. ex- 
ponential case (\ = 5 — 0.2 ). The "■ ■ ■ " (dashed line) indicates true value of the parameter 
in each case. 



the initial parameter estimates for the distribution of X and T using the 
right truncated distribution. These initial estimates are way off in all cases, 
which explains why standard practice is unsatisfactory in this nontrivial sit- 
uation. We summarize our simulation result in Table 1. The first two rows 
in Table 1 are of special interest since we assumed X,T ~' Exp(A = 6). 
As shown in Appendix B.5, the exact likelihood given in equation (2.1) can 
be solved numerically in this case. For Tq = 6 the exact likelihood based 
MLE yields A = 0.22 with asymptotic standard deviation ax = 0.026. For 
To = 5, we get A = 0.18 with asymptotic standard deviation a\ = 0.028. In 
both of these cases our simulation result is very close to the true value 
{X = 5 = 0.2) even though we did not use the information that X = 6 in our 
proposed algorithm. In Figure 4 we present pictorially the result for these 
two cases. This supports the viability of our algorithm. Next we explore 
non-i.i.d. cases. Figure 3 presents the case for X Exp(A) and T ~ Exp((5) 
with two different observation times (Tq = 4, 6). In the first case, we choose 
the true model parameters in such a way that about 50% of the cases are 
observed (i.e., C > 100). Figure 3(a) and (b) present the case when A = 0.5 
and 6 = 0.2. We observe 108 and 66 units for Tq = 6 and 4, respectively. As 
expected, the case with more units produces better estimates. Nevertheless, 
we point out that for Tq = 4, even though we observe only about 33% of 
the units, the final parameter estimates are still noticeably close to the true 
parameter values. Similar observations could be made for the other choice 
of parameter values in Figure 3(c) and (d). To elucidate the problem when 
using the exact maximum likelihood based approach, we have also plotted 
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(c) (d) 

Fig. 5. Log-likelihood surface plots for the Exponential-Exponential model with A = 0.5 
and 5 = 0.2. (a,) is obtained when we use all the observations (C = 2Q0). (h) is obtained 
when To = 6 and C = 108. (c) is obtained when To = 4 and C = 66. Likelihood becomes 
flatter as C ^, thus making MLE search a difficult task. (A) is obtained for a specific 
iteration when imputation is used (97 imputed samples) for To = 4 and C = 66. 

the log-likelihood surface (obtained via equation (2.1) and numerical inte- 
gration) in Figure 5 for the case A = 0.5 and 5 = 0.2. Figure 5(a) represents 
the case when we have complete observations for all units (C = A^). How- 
ever, as To shrinks, C goes down, and, as a result, the likelihood surface 
becomes very flat. Hence, searching for the MLE becomes computationally 
challenging and often leads to large variance. We have noted this problem 
earlier in Section 2. Figure 5(d) presents the log- likelihood surface obtained 
via equation (6.1) when imputation is in use. This representative plot is ob- 
tained for a specific iteration when 97 units are imputed while running the 
algorithm described in Section 5. The flatness of the resulting log-likelihood 
surfaces in Figure 5(c) and (d) is an indicator of computational difficulties 
in finding the MLE for each case. Next, in Figure 6 we describe the iteration 
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X estimate plot 



T (shape) estimate plot 




200 400 600 800 1000 200 400 600 800 1000 

Iteration Iteration 



(a) 



(b) 



T (scale) estimate plot 




1^ 




200 400 600 800 1000 

Iteration 



(c) 

Fig. 6. Plot of the maximum likelihood estimate for the Exponential- Weihull model over 
different iterations. The (dashed line) indicates true value of the parameter in each 



result when X ~ Exp(A) and T ~ Weibull(/3, 6). In Figure 7 we describe the 
iteration result when X ~ Weibull(/3, 0) and T~Exp(A). In all cases the 
final estimates are quite close to the true model parameters. Though not 
reported here, we obtain similar results with the gamma distribution. For 
details of the sampling from a truncated gamma distribution, please refer 
to Damien and Walker (2001). We have confined our simulation exploration 
only to commonly used reliability distributions; however, we are hopeful 
that the algorithm presented here will also work for other distributions with 
nonnegative support. 
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200 400 600 800 1000 

Iteration 



(c) 

Fig. 7. Plot of the maximum likelihood estimate for the Weibull-Exponential model over 
different iterations. The (dashed line) indicates true value of the parameter in each 

case. 

8. Motivating application. The data set that we will analyze using the 
current procedure came from an industrial house producing residential fur- 
nace components during one week in May 2001. We consider a batch with 
N = 400 units. The data consist of C = 133 pairs of points as observed units 
(i.e. , {xi,ti ) , which have failed within the observation time of seven years 
from the date of manufacturing. Figure 8(a) shows a violin plot for instal- 
lation and failure times. The violin plot is a combination of a box plot and 
a kernel density plot. There is no specific information available about the 
remaining units. We are assuming that there exists no unit which has failed 
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X estimate plot 



Parallel Violin plots for the Installation and Failure time 





(a) 

T (shape) estimate plot 



(b) 

T (scale) estimate plot 



400 600 
Iteration 




(d) 



Fig. 8. Plot (a.) represents violin plot for the observed 133 units. Other plots represent 
the maximum likelihood estimates for the Furnace data over different iterations. The "■ ■ ■ " 
(dashed line) indicates finally estimated mean of the parameter in each case. 



but was not reported. In practice, this could have happened for many other 
reasons. In the present context the rehabihty engineers beheve that it is ap- 
propriate to model installation time {X) using an exponential distribution, 
while failure time (T) is modeled according to a Weibull distribution [Jager 
and Bertsche (2004); Zhu (2007)]. It should be noted that seasonality plays 
an important role in selling, installation and duty cycles (how rigorously 
the unit is being used) of the product. However, since in the present case 
we consider only a single batch, we assume that these effects will be simi- 
lar for every unit in the batch. When comparing the units produced under 
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different batches (and possibly produced at different times of the year), ad- 
ditional care is required as the independence assumption between X and T 
becomes questionable. This is due to the fact that some installation times 
are associated with severe duty cycles and more reliability problems. 

Before running the algorithm we divide the installation times as well as 
failure times by their corresponding standard deviation estimated from 133 
samples. This rescaling is done for numerical stabilization only, which re- 
sults in faster convergence of the algorithm. Rescaled random variables have 
straightforward relationships with the original variables, without any drastic 
change to the distributional form. We run the algorithm for 1000 iterations, 
however, convergence (with p = 5, e = 0.0005) was achieved much earlier. 
We discard the first 100 iterations as burn-in and report the estimates on 
the basis of the remaining 900 iterations in Table 2. For model comparison 
purposes we have also investigated separately the case where T is assumed 
to follow the exponential distribution, without altering the distribution of 
X. In each case we obtain the initial parameter estimates using the right 
truncated distributions. Figure 8 represents the case for the Exponential- 
Weibull model combination. Though the Exponential-Exponential model 
parameter is different from the previous choice (see Table 2), the density 
plot of the two distributions of T are quite similar as depicted in Figure 
9(b). We have also compared the predictive performance of different mod- 
els in Figure 9(c), including the usual practice of truncated distributions 
without any imputation. We estimated the expected number of failures to 
be observed for different observation times over an interval of six months. 
This expected failure number is then compared with the observed failure 
number for the current data set. This required repeated re-estimation of 
model parameters at different time points. As can be seen, the truncated 
models have a huge overestimation problem throughout the study period. 
This again justifies our earlier criticism of current practice. Imputed models 
produce stable estimates and do much better even at the very early stage 
of product lifetime with only limited data. The Exponential-Weibull model 
choice does a little better than the Exponential-Exponential model. How- 
ever, they are very much comparable as expected from Figure 9(b). It is 
desirable to estimate the expected failure number accurately for two main 
reasons. First, by accurately estimating warranty claims, an estimate of re- 
quired financial reserves can be performed. This has immense implications in 
terms of future financial resource management. Second, it is desired to con- 
tinuously improve the quality of consumer products, especially at the very 
high quality levels enjoyed by many consumer products today. All these as- 
pects necessarily depend upon the accurate and efficient estimation of the 
reliability parameters (in X and T). The method described in this paper 
provides a first step in this direction. 
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Table 2 

Estimates for N = 400 units m a single batch. C — 133 units have complete observations 





Initial 




Simulation 


Average No. 


Convergence 


Time in 


Distribution 


estimate 




result 


imputations 


iteration 


second 


X ~ Exp(A) 


A = 0.9 


A = 


= 0.57, CTA = 0.014 


260 


167 


381 


T ~ Weibull(/3, e) 


/3 = 0.6, 




= 0.81, ?^ =0.004 










61 = 3.18 


e-- 


= 14.47, ?e = 0.4 








T ~ Exp(5) 


5 = 0.51 


5 = 


0.079, ?a =0.001 


263 


45 


421 



9. Concluding remarks. Unlike electronic commodities, item specific track- 
ing is not a feasible solution for many large scale industrial operations. 
Hence, the availability of both "complete" and "partial" information is quite 
common. In addition, except for very rare occasions, there are hardly any 
situations where all units in a batch start working at the same time. Un- 
availability of the installation time in a timely fashion is a major challenge 
to reliability engineers. Because of confidentiality issues we can not reveal 
any company specific information. However, we would like to mention that 
the above problem exists in different industrial sectors, and there is no clear 
solution thus far. In this paper we have proposed a computational approach 
to solve the problem with the optimal usage of partial and complete infor- 
mation. From a reliability engineer's perspective, this current approach is 
simple, fast and also has straightforward interpretability. 

The primary focus of any reliability analysis is the failure time. How- 
ever, the waiting time for the installation is also very important in the sense 
that it provides valuable market specific information from the sales per- 
spective, including seasonality and periodic sales patterns. In our approach 
we have targeted simultaneous estimation for both installation and failure 
time parameters in a combined fashion. To the best of our knowledge, this 
is the first attempt to do so. Finally, we would like to point out some of 
the assumptions that we have made in this paper, a violation of which will 
require more research. First, we have assumed that installation time and 
failure time are independent. This may be questionable in some situations 
as discussed in Section 8. Second, there is no aging effect for the units in- 
stalled at different time points. Finally, we made the assumption that the 
distributional form of both installation and failure times is known. While 
for most of the legacy industrial products, in-house experts have a good idea 
about this from historical knowledge, it is of theoretical interest to see the 
effect of convergence and the quality of parameter estimates under incorrect 
parametric model specification. One way to avoid this is to choose a larger 
class of models. From the reliability perspective there is considerable effort 
to generalize Weibull and other popular reliability distributions [see Bali 
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T estimate plot Distribution of the scaled faiiure time 




(a) (b) 



Performance of different model 



Model choice 
Observed Failure 
Expo-Weibull {with imputation) 
Expo-Expo (with imputation) 
Expo-Weibull (truncated) 
Expo-Expo (truncated) 



Months ahead (in year) 

(c) 



Fig. 9. On the left, (&), maximum likelihood estimate for the furnace data when 
T ~ Exponential distribution. In the middle, (h), density plot for two different model 
choices for T . Both look similar. On the right, (c), it represents performance of different 
models compared with the observed failure. Truncated distribution with no imputation per- 
forms very poorly with huge overestimation. Performance of imputed models are far better 
and the Exponential- Weibull model choice does the best job. 



(2003) and Shao (2004)]. However, the resultant estimation procedure will 
be more involved. Another possibility is a nonparametric extension; however, 
the resulting procedure will be much more complex. In an ongoing work we 
are also exploring the exact probabilistic and inferential procedure based on 
equation (2.1). 
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APPENDIX A: PROOF OF THEOREM 3.2. 
We can use the inequality (3.4) to argue that the following holds: 

St{To - Xk) < P[T > To - X\xk <X< xk+i] < St{To - Xk+i), 
St{To - Xk-i) < P[T > To - X\xk-i <X<Xk]< St{To - Xk). 

Combining both of these yields the proof. 

APPENDIX B: MAXIMUM LIKELIHOOD ESTIMATION. 

We concentrate here on Exponential and Weibull distribution as used 
in the simulation, though other distributions with positive support, such 
as gamma and log-normal, can also be considered. Most of the results are 
published elsewhere and referenced as required. 

B.l. Truncated exponential. Let X ~ Exp(A) with < X < Tq. The 
p.d.f. is given by 

ff IX ^ N Aexp(-xA) 



1 - exp(-ToA) 



If we have n observations, then differentiating the log-likelihood equation 
with respect to A and equating it to zero yields 

1 roexp(-roA) _ 

7 — ;^rTT — x = U. 

A 1 - exp (-TqA) 

The above equation needs to be solved numerically to get the MLE of A. 

B.2. Randomly right censored exponential. Let T ~ Exp (A) and we ob- 
serve T* = min{r, Cj.}, where in the current context Cr = Tq — X and X is 
another random variable denoting installation time. Let us denote our sam- 
ples as {t* ,5i}^^i, where 6i = 1 means the sample is an actual observation 
and means it is censored. If we have Y17=i ^i — C true observations, then 
the log-likelihood is given by 

C n-C 

L{X) = clogX - X^ti - Xj2iTo - xj), 
i=i j=i 



which upon equating to yields A : 



c 
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B.3. Truncated Weibull. The MLE calculation for the truncated Weibull 
distribution is somewhat involved and may not always exist. Some explicit 
mathematical formulations with the required regularity conditions are de- 
scribed in Mittal and Dahiya (1989). We briefly mention only the final result 
here that has been used in this paper. Suppose X ~ Weibull(/3, but with 
< X < To. Let us denote by Y = ^. Unfortunately, the MLE for /3 is 
not available in closed form and needs to be solved numerically using the 
equation 



n 



^//5+Er=iiog?/i 



+ 



exp 



1 



0. 



Once we know /?, the MLE of 9 is 



e = Tn 



B.4. Randomly right censored Weibull. Suppose T ~ Weibull(/3, 9) . Sim- 
ilar to the randomly right censored exponential case T* = minjT, Cr}, where 
in the current context Cj- = Tq — X . We denote our data set as {t*,5i}f^i 
and EiLi — The MLE is given explicitly in Shao (2004) and Lemon 
(1975), which again needs to be solved numerically for (3 using the equation 



1 I ELi'^aogt: 

/? c 



EiLi(^n^iog^: 
ELi(in^ 



0. 



Once we know /3, the MLE of 9 is 

c 



B.5. Derivation of the exact MLE for i.i.d. exponential case. We assume 



,i.i.d 



X,T '~'Exp(A). The complete likelihood is given by 
L(A) oc (A)2^e-^S?=i(^»+*')[e~^^« +roAe- 



XToiN-C 



Now differentiating the log-likelihood equation with respect to A and equat- 
ing it to zero yields 



The equation needs to be solved numerically for A to obtain MLE. 
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SUPPLEMENTARY MATERIAL 

Furnace Data Set and R Code for Furnace Data as well as Simulation 
for all Models Considered in the Paper (DOl: 10. 1214/10- AOAS348SUPP; 
.zip). R code is used for the simulation as well as real data analysis. 

Supplementary material has five files: 

1. Furnace data in MS Excel format (data.xls). 

2. Code for analyzing furnace data (code_furn.doc). 

3. Code for the Exponential-Exponential model (new_code_Exp(2).doc). 

4. Code for the Exponential- Weibull model (new_code_ExpWeb.doc). 

5. Code for the Weibull-Exponential model (new_code_WebExp.doc). 

For the simulation examples data sets are generated on the fly at the begin- 
ning of the code. No special R package is required to run the codes. All the 
codes are commented for the ease of understanding. 
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